// This code converts AlpGen's mode 2 unw files to a lhe file. 

// The produced lhe file can be converted into a root file using LHEFConvert 
// included into the ExRoot package.
// This version comlpetely fakes color and mother-child informations. 
// On the contrary, weights and scales information are copied from the original unw.
// Future version will properly deal with these data
// currently it is just good for kinematic studies at partonic level with no 
// showering or hadronization


#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iterator>
#include <vector>
#include <iostream>
#include <sstream>
#include <string>
#include <cstring>
#include <cstdlib>
#include "math.h"
#include <stdio.h>

using namespace std;

int main(int argc, char* argv[]) {

	if (argc < 3) {
		cout << " Usage: unw2lhe input_file1 input_file2 ... input_fileN output_file" << endl;
		cout << " input_file - input file in UNW format," << endl;
		cout << " output_file - output file will be in LHE format." << endl;
		return 1;
	}
	for (int ic = 0; ic < argc - 1; ic++) {
		if (strcmp(argv[argc - 1], argv[ic])==0) {
			cout << "inpunt and outfile cannot be the same" << endl;
			return 1;
		}
	}
	int ninout = 0;
	vector<string> lines;
	vector<string> fields;
	string eventweight;
	string pdfscale;
	string processID;

	ofstream lhe;
	lhe.open(argv[argc - 1]);
	string opening;
	opening = "<LesHouchesEvents version=\"1.0\">\n<!--\n-->\n<init>\n  2212  2212  0.70000000000E+04  0.70000000000E+04 0 0 10042 10042 3 1\n  0.11750367120E-02  0.10499153982E-04  0.23730090000E-05 661\n</init>";
	cout << "The resulting LHE file will have the following Opening" << endl;
	cout << opening << endl;
	lhe << opening << endl;
	lhe.close();
	string temp;
	string tempfield;
	istringstream infield;
	istringstream instream;
	int blockposition = 0;
	string in_pid;
	float in_p_L;
	float in_E;
	string pid;
	float pX;
	float pY;
	float pL;
	float m;
	float E;

	// LOOP ON FILES
	for (int infilecounter = 1; infilecounter < argc - 1; infilecounter++) {
		ifstream unw(argv[infilecounter]);
		cout << "Processing " << argv[infilecounter];

		//LOOP  on LINES
		while (getline(unw, temp)) {
			//cout<<temp<<endl;
			lines.push_back(temp);
			instream.clear();
			instream.str(temp);
			fields.clear();

			if (blockposition == 0) {

				if (ninout == 0) {
					while (getline(instream, tempfield, ' ')) {

						if (tempfield != "" && tempfield != " ") {
							fields.push_back(tempfield);
							//cout<<tempfield<<endl;
						}
						//cout<<fields[0]<<endl;
					}
					//all fields are now filled
					// read the number of particles
					// and the event informations written in the same line
					processID = fields[1];
					ninout = atoi(fields[2].c_str());
					eventweight = fields[3];
					pdfscale = fields[4];

				}

				// event first line
				// n. of particles,  process id, weight, scale, aEM, Astrong
				lhe.open(argv[argc - 1], ios::app);
				lhe << "<event>\n" << ninout << " " << processID << " " << eventweight << " " << pdfscale << " 0.7816531E-02  0.1300000E+00" << endl;
				lhe.close();

			}
			//incoming particle lines
			if (blockposition <= 2 && blockposition != 0) {

				while (getline(instream, tempfield, ' ')) {

					if (tempfield != "" && tempfield != " ") {
						fields.push_back(tempfield);
						//cout<<tempfield<<endl;
					}
					//cout<<fields[0]<<endl;
				}
				//all fields are now filled

				in_pid = fields[0];
				//cout<<in_pid<<endl;
				in_p_L = atof(fields[3].c_str());
				in_E = sqrt(in_p_L * in_p_L);

				lhe.open(argv[argc - 1], ios::app);
				lhe << "       " << in_pid << "   -1    0    0  502  501  0.00000000000E+00  0.00000000000E+00  " << in_p_L << "  " << in_E << "  0.00000000000E+00 0. -1." << endl;
				lhe.close();
				fields.clear();

			}
			//outgoing particles
			if (blockposition > 2) {
				while (getline(instream, tempfield, ' ')) {
					if (tempfield != "" && tempfield != " ") {
						fields.push_back(tempfield);
					}
					//fields are ready to be processed and outputted to the lhe file
				}

				pid = fields[0];
				pX = atof(fields[3].c_str());
				pY = atof(fields[4].c_str());
				pL = atof(fields[5].c_str());
				m = atof(fields[6].c_str());
				E = sqrt(pX * pX + pY * pY + pL * pL + m * m);
				lhe.open(argv[argc - 1], ios::app);
				lhe << "       " << pid << "    1   11   11    0  505 " << pX << " " << pY << " " << pL << " " << E << "  " << m << " 0.  1." << endl;
				lhe.close();
				fields.clear();

			}
			blockposition += 1;
			if (blockposition == ninout + 1) {
				lhe.open(argv[argc - 1], ios::app);
				lhe << "</event>\n" << endl;
				lhe.close();
				blockposition = 0;
				ninout=0;
			}
			/*
			 ostream_iterator<string> os_it(cout," ");
			 copy(fields.begin(),fields.end(),os_it);
			 cout<<"-----------line has been processed"<<endl;
			 */
		}
	}
	lhe.open(argv[argc - 1], ios::app);
	lhe << "</LesHouchesEvents>\n" << endl;
	lhe.close();
	return 0;
}
